/* yyt05/cpp
 * Checking of Yuan, Yuan, Ting research from PRB 71, 104505 (2005)
*/

/* STANDARD LIBRARY */
#include<cstdio>
#include<stdlib.h>
#include<iostream>
#include<string>
#include<math.h>

/* ADDITIONAL LIBRARY */
#include"constants.h" // constants
#include"supToGSL.h" // support to GSL library
#include"yyt05Funcs.h" // supprot to main program (additional functions)
#include"yyt05Set.h" // set of functions to solve using GSL
#include"colorText.h" // functions which change color of text (for fun)
#include"jkClass.h"
#include"thread.h"
using namespace std;

int yyt2dGA (const gsl_vector *_x, void *_Par, gsl_vector * fun)
{
	Parameters *Par = reinterpret_cast<Parameters*>(_Par);

	const unsigned int N = 5; // we have set of 5 functions to solve
	double x[N]; // vector of arguments
	for(unsigned int i=0; i<N; i++) { x[i] = gsl_vector_get(_x, i); }
	Parameters *X0 = new Parameters(N, x);
	yytGA *yyt = new yytGA(Par, N, X0);

	double *y = yyt->sumationAll();

	printf("--->x = ");
	for(unsigned int i=0; i<N; i++) printf("%20.16f ", x[i]);
	printf("\n--->y = ");
	for(unsigned int i=0; i<N; i++) printf("%20.16f ", y[i]);
	printf("\n\n");

	// I put finally result from y to gsl vector "fun"
   for(unsigned int i=0; i<N; i++) { gsl_vector_set (fun, i, y[i]); }
	// cerr << "DEBUG in yyt2dGA yyt " << yyt  << endl;
   delete yyt;
   delete X0;
	// cerr << "DEBUG in yyt2dGA" << endl;
   return GSL_SUCCESS;
}

void* Run(void *wsk)
{
	Clock *clock = new Clock();
	ParThre *parThre = static_cast<ParThre*>(wsk);
	int id = parThre->id;
	double *x0 = parThre->x0;
	double *par = parThre->par;
	unsigned int nPar = parThre->nPar;
	unsigned int d = parThre->d;
	unsigned int nRepetitions = parThre->nRepetitions;
	double start = parThre->start;
	double end = parThre->end;
	const char *file = static_cast<const char*>(parThre->file);
	const char *title = static_cast<const char*>(parThre->title);
	const size_t nEquations = parThre->nEquations;
	const unsigned int MaxNumberOfShake = parThre->MaxNumberOfShake;
	unsigned maxInterp = parThre->maxInterp;
	int (*func)(const gsl_vector*, void*, gsl_vector*) = parThre->func;
	pthread_mutex_t *mut = parThre->mut;

	FILE *fp = fopen(file, "w");
	if(fp==NULL)
	{
		fprintf(stderr, "Error with opening %s to write!\n", file);
		exit(EXIT_FAILURE);
	}

	Parameters Par(nPar, par); // list of parameters // 7 or 9!!!
	Parameters X0(nEquations, x0); // start point
	MultirootSolverF S(func, nEquations); 	// which function will be solved

	fprintf(fp, "# sdelta=%.3f, J=%.3f, t=%.3f, U=%.3f", par[0], par[1], par[2], par[3]);
	fprintf(fp, " beta=%.3f, L=%d\n", par[4], (int)par[5]);
	fprintf(fp, "%s", title);

	double **y = new double*[nEquations];
	for(unsigned int i=0; i<nEquations; i++) { y[i] = new double[nRepetitions+1]; }
	Parameters *Result = NULL;
	unsigned int nShake = 0;
	bool isStuck = false;
	double *x = new double[nRepetitions+1];

	for(unsigned int i=0; i<=nRepetitions; i++)
	{
		// a critical section required - start
		mutexClose(mut);
//		cleanLine(id+1);
//		gotoxy(0,id+1);
		float bar = static_cast<float>(i)/static_cast<float>(nRepetitions);
		printf("|");
		for(unsigned int j=0; j<50; j++)
		{
			if(bar*50. <= static_cast<float>(j)) printf(" ");
			else printf("=");
		}
		printf("| ");
		int iswitch = i%4;
		switch(iswitch)
		{
			case 0: printf("|"); break;
			case 1: printf("/"); break;
			case 2: printf("-"); break;
			case 3: printf("\\"); break;
		}
		if(isStuck)
		{
			printf(" %.0f%%", bar*100.);
			printf("  aborted :-(\n");
			mutexOpenUnlock(mut);
			for(unsigned int i=0; i<nEquations; i++) { delete []y[i]; }
			delete []y;
			delete clock;
			fclose(fp);
			return NULL; }
		if(bar < 1)
		{
			printf(" %.0f%%", bar*100.);
			float time = clock->estimateEnd(bar);
			const char *descriptorT = (const char*)"sec";
			if(time>60.) { time /= 60.; descriptorT = (const char*)"min"; }
			else if(time>3600) { time /= 3600.; descriptorT = (const char*)"hours"; }
			printf(" end in %.0f %s\n", time, descriptorT);
		}
		else { printf("  complete :-)\n"); }
		mutexOpenUnlock(mut);
		// a critical section required - end

		double h = (end-start)/static_cast<double>(nRepetitions);
		if(nRepetitions==0) h=0.;
		Par[d] = static_cast<double>(i)*h + start;
		x[i] = Par[d];
		if(DEBUG) X0.print();
		// cerr << "DEBUG1" << endl;
		S.solveFor(&Par, &X0);
		// cerr << "DEBUG2" << endl;
		Result = S.getResult();

		if(Result == NULL) // shake of system
		{
			fprintf(stderr, "%d: NULL error for %.6f (%d time)\n", id, Par[d], nShake);
			X0.shake(++nShake);

			if(nShake>MaxNumberOfShake)
			{
				fprintf(stderr, "error for: Par[%d] = %.6f\n", d, Par[d]);
				fprintf(stderr, "It is imposible to solve the set of equations, sorry\n");
				isStuck = true;
			}

			i--;
			continue;
		}

		for(unsigned int j=0; j<nEquations; j++) { y[j][i] = (*Result)[j]; }

//		for(unsigned int j=0; j<nEquations; j++) { X0[j] = (*Result)[j]; }
		X0.interpolation(x, y, maxInterp, i, x[i], h);

		fprintf(fp, "%.9f\t", Par[d]);
		S.printResult(fp);
		nShake = 0;
		Result->deleteP(); delete Result;
	}

	for(unsigned int i=0; i<nEquations; i++) { delete []y[i]; }
	delete []y;
	delete clock;
	fclose(fp);
	return NULL;
}

int main()
{
	unsigned int nThreads=1; // number of threads
	pthread_t *id[nThreads];
	pthread_mutex_t *mut = mutexNew();
	ParThre **parThre = new ParThre*[nThreads];
	parThre[0] = new ParThre();

	// XV.3 2D, standard, testy rozdzielczosci: wykres od L;
	parThre[0]->MaxNumberOfShake = 2; // how many we can shake system?
	parThre[0]->maxInterp = 3; // how many points will be use to interpolation
	parThre[0]->mut = mut;
	parThre[0]->d = 3; // which parameter will be changed
	parThre[0]->nRepetitions = 100; // and how many times
	parThre[0]->start = 0.; // from which value
	parThre[0]->end = 20; // to which value
	// SGA
	parThre[0]->nEquations = 5; // how many equations
	parThre[0]->func = &yyt2dGA;
	parThre[0]->title = "# 2D: t-J-U: GA - wykres od sdelta\n\
			# U\tstatus\titer\tBDelta\t\tchi\t\tm\t\td\t\tmu\n";
	const unsigned int nPar = 6; // 6 or 8
	const unsigned int nVar = 5;
	parThre[0]->nPar=nPar;
	double par[nThreads][nPar];
	double x[nThreads][nVar];

	unsigned int i=0;
	// 0
	parThre[i]->id = i; // id of current thread
	par[i] = {0.0, 0.333, 1., 8., 500., 512.};
	x[i] = {0.004, 0.40, 0.012, 0.17, -0.19};
	parThre[i]->par = par[i]; // parameters: sdelta, J, t, U, beta, L
	parThre[i]->x0 = x[i]; // start point: // BDelta, chi, m, d, mu
	parThre[i]->file = "plots/tJU/tjuGA_sdelta0_U0-20.dat";
	id[i] = threadCreate(Run, parThre[i]);


	// Running
	mutexClose(mut);
//	cleanScrean();
	textRed("Program MROS - test t-J-U dla 2 wymiarow, GA\n");
	mutexOpenUnlock(mut);

	for(unsigned int i=0; i<nThreads; i++)
	{
		//if(*id[i] == NULL) { continue; }
		fprintf(stderr, "[[%d:", i);
		threadJoinIgn(*id[i]);
		fprintf(stderr, "%d]]", i);
	}

	// Cleaning
//	gotoxy(0,nThreads+2);
	mutexDestroy(mut);
	for(unsigned int i=0; i<nThreads; i++) { delete parThre[i]; }
	delete[] parThre;
	return 0;
}
//:~ end of main
